﻿using System;
using System.Globalization;

namespace Kyoh
{
//	[Obsolete("未完成なため使用できません。", true)]
	/// <summary>
	/// 一般行列を扱うクラスです。
	/// </summary>
	public class Matrix16 : ICloneable
	{
		#region フィールド

		/// <summary>
		/// 行列の内部表現である二次配列を保持するフィールドです。
		/// </summary>
		int[,]	
					mValue,
					mTempMat1,
					mInverse,
					mMat_q,
					mMat_r,
					mMat_g;

		int
					mRow,
					mCol,
					mLowerSize;

		public const int PointPosition = 16;

		#endregion

		#region プロパティ・インデクサ

		/// <summary>
		/// インスタンスが持つ行列の行数を取得します。
		/// </summary>
		public int Row
		{
			get { return mRow; }
		}

		/// <summary>
		/// インスタンスが持つ行列の列数を取得します。
		/// </summary>
		public int Col
		{
			get { return mCol; }
		}

		/// <summary>
		/// 指定された行・列の位置にある値を取得・設定します。
		/// </summary>
		/// <param basename="row">行位置を指定します。</param>
		/// <param basename="col">列位置を指定します。</param>
		/// <returns>指定された行・列の位置にある値を返します。</returns>
		public double this[int row, int col]
		{
			get { return mValue[row, col] / (double)(1<<PointPosition); }
			set { mValue[row, col] = (int)(value * (1<<PointPosition)); }
		}

		public int GetValueAt(int row, int col)
		{
			return mValue[row, col];
		}

		public void SetValueAt(int row, int col, int value)
		{
			mValue[row, col] = value;
		}

		#endregion

		#region コンストラクタ

		/// <summary>
		/// 行列クラスのインスタンスをデフォルトで生成します。
		/// このとき、生成されたインスタンスは1*1の行列になります。
		/// </summary>
		public Matrix16() : this(1, 1) { }

		/// <summary>
		/// 指定された数の次元で正方行列のインスタンスを生成します。
		/// このとき、生成されたインスタンスは<paramref basename="total"/>*<paramref basename="total"/>の行列になります。
		/// </summary>
		/// <param basename="total">次元数を指定します。</param>
		public Matrix16(int size) : this(size, size) { }

		/// <summary>
		/// 指定された行数・列数で行列のインスタンスを生成します。
		/// このとき、生成されたインスタンスは<paramref basename="row"/>*<paramref basename="col"/>の行列になります。
		/// </summary>
		/// <param basename="row">行数を指定します。</param>
		/// <param basename="col">列数を指定します。</param>
		public Matrix16(int row, int col)
		{
			mRow = row;
			mCol = col;
			if (mRow < mCol)
			{
				mLowerSize = mRow;
			}
			else
			{
				mLowerSize = col;
			}
			mValue = new int[row, col];
			mMat_q = new int[mLowerSize, mLowerSize];
			mMat_r = new int[mLowerSize, mLowerSize];
			mMat_g = new int[mLowerSize, mLowerSize];
			mInverse = new int[col, row];
			mTempMat1 = new int[mLowerSize, mLowerSize];
		}

		/// <summary>
		/// 指定された倍精度浮動小数点数型二次配列を内部表現に持つ行列のインスタンスを生成します。
		/// このとき、生成されたインスタンスはmatrix[a,b]に対してa*bの行列になります。
		/// </summary>
		/// <param basename="matrix">内部表現となる倍精度浮動小数点数型二次配列を指定します。</param>
		public Matrix16(double[,] matrix)
			: this(matrix.GetLength(0), matrix.GetLength(1))
		{
			int i, j;
			for (i = 0; i < mRow; i++)
				for (j = 0; j < mCol; j++)
					mValue[i, j] = (int)(matrix[i, j] * (1 << PointPosition));
		}

		/// <summary>
		/// 指定された小数点位置<see cref="Matrix16.PointPos"/>の固定小数点数型二次配列を内部表現に持つ行列のインスタンスを生成します。
		/// このとき、生成されたインスタンスはmatrix[a,b]に対してa*bの行列になります。
		/// </summary>
		/// <param basename="matrix">内部表現となる倍精度浮動小数点数型二次配列を指定します。</param>
		public Matrix16(int[,] matrix)
			: this(matrix.GetLength(0), matrix.GetLength(1))
		{
			int i, j;
			for (i = 0; i < mRow; i++)
				for (j = 0; j < mCol; j++)
					mValue[i, j] = matrix[i, j];
		}

		#endregion

		#region 単項演算子用
		/// <summary>
		/// このインスタンスを正負反転した行列のインスタンスを取得します。
		/// </summary>
		/// <returns>このインスタンスを正負反転した行列のインスタンスを返します。</returns>
		public Matrix16 Negative()
		{
			double[,] ret = new double[mRow, mCol];
			int i, j;
			for (i = 0; i < mRow; i++)
				for (j = 0; j < mCol; j++)
					ret[i, j] = -mValue[i, j];
			return new Matrix16(ret);
		}
		#endregion

		#region 四則演算
		/// <summary>
		/// 行列の加算演算を定義します。
		/// ここでの加算演算によりこのインスタンスが変化することはありません。
		/// 加算結果が新しいインスタンスに格納されて返ります。
		/// </summary>
		/// <param basename="mat">このインスタンスに足す行列のインスタンスを指定します。</param>
		/// <returns>加算演算の結果を返します。</returns>
		public Matrix16 Add(Matrix16 mat)
		{
			if (mRow != mat.Row || mCol != mat.Col)
				return null;
			double[,] ret = new double[mRow, mCol];
			int i, j;
			for (i = 0; i < mRow; i++)
				for (j = 0; j < mCol; j++)
					ret[i, j] = mValue[i, j] + mat.mValue[i, j];
			return new Matrix16(ret);
		}
		/// <summary>
		/// 行列の減算演算を定義します。
		/// ここでの減算演算によりこのインスタンスが変化することはありません。
		/// 減算結果が新しいインスタンスに格納されて返ります。
		/// </summary>
		/// <param basename="mat">このインスタンスから引く行列のインスタンスを指定します。</param>
		/// <returns>減算演算の結果を返します。</returns>
		public Matrix16 Sub(Matrix16 mat)
		{
			if (mRow != mat.Row || mCol != mat.Col)
				return null;
			double[,] ret = new double[mRow, mCol];
			int i, j;
			for (i = 0; i < mRow; i++)
				for (j = 0; j < mCol; j++)
					ret[i, j] = mValue[i, j] - mat.mValue[i, j];
			return new Matrix16(ret);
		}
		/// <summary>
		/// 行列の乗算演算を定義します。
		/// ここでの乗算演算によりこのインスタンスが変化することはありません。
		/// 乗算結果が新しいインスタンスに格納されて返ります。
		/// </summary>
		/// <param basename="mat">このインスタンスに乗ずる行列のインスタンスを指定します。</param>
		/// <returns>乗算演算の結果を返します。</returns>
		public Matrix16 Mult(Matrix16 mat)
		{
			if (mCol != mat.mRow)
				return null;
			double[,] ret = new double[mRow, mat.mCol];
			int i, j, k;
			for (i = 0; i < mRow; i++)
				for (j = 0; j < mat.mCol; j++)
					for (k = 0; k < mCol; k++)
						ret[i, j] += mValue[i, k] * mat.mValue[k, j];
			return new Matrix16(ret);
		}
		/// <summary>
		/// 行列の除算演算を定義します。
		/// ここでの除算演算によりこのインスタンスが変化することはありません。
		/// 除算結果が新しいインスタンスに格納されて返ります。
		/// </summary>
		/// <param basename="mat">このインスタンスから除す行列のインスタンスを指定します。</param>
		/// <returns>除算演算の結果を返します。</returns>
		public Matrix16 Div(Matrix16 mat)
		{
			return Mult(mat.Inverse());
		}
		#endregion

		#region 転置行列
		/// <summary>
		/// このインスタンスの行列の転置行列を取得します。
		/// </summary>
		/// <returns>インスタンスの転置行列を返します。</returns>
		public Matrix16 Transposed()
		{
			double[,] ret = new double[mCol, mRow];
			int i, j;
			for (i = 0; i < mRow; i++)
				for (j = 0; j < mCol; j++)
					ret[j, i] = mValue[i, j];
			return new Matrix16(ret);
		}
		/// <summary>
		/// 指定されたインスタンスの行列の転置行列を取得します。
		/// </summary>
		/// <param basename="sqm">転置行列を得たい行列のインスタンスを指定します。</param>
		/// <returns>指定されたインスタンスの行列の転置行列を返します。</returns>
		static public Matrix16 TransposedOf(Matrix16 sqm)
		{
			double[,] ret = new double[sqm.Col, sqm.Row];
			int i, j;
			for (i = 0; i < sqm.Row; i++)
				for (j = 0; j < sqm.Col; j++)
					ret[j, i] = sqm[i, j];
			return new Matrix16(ret);
		}
		#endregion

		#region 逆行列

		//UNDONE: リファクタ

		/// <summary>
		/// 一般逆行列を算出します。
		/// </summary>
		/// <returns>(TG * G)~ * TG</returns>
		public Matrix16 Inverse()
		{
			int i,
				j,
				k;

			// 転置行列(左)の本行列(右)の乗算行列を求めます
			if (mRow < mCol)
			{
				for (i = 0; i < mRow; i++)
				{
					for (j = 0; j < mRow; j++)
					{
						mTempMat1[i, j] = 0;
						for (k = 0; k < mCol; k++)
							mTempMat1[i, j] += mValue[i, k] * mValue[j, k];
					}
				}
			}
			else if (mCol < mRow)
			{
				for (i = 0; i < mCol; i++)
				{
					for (j = 0; j < mCol; j++)
					{
						mTempMat1[i, j] = 0;
						for (k = 0; k < mRow; k++)
							mTempMat1[i, j] += mValue[k, j] * mValue[k, i];
					}
				}
			}
			else
				for (i = 0; i < mRow; i++)
					for (j = 0; j < mCol; j++)
						mTempMat1[i, j] = mValue[i, j];

			#region 正則行列の逆行列化

			for (i = 0; i < mLowerSize; i++)
			{
				for (j = 0; j < mLowerSize; j++)
					mMat_g[i, j] = mTempMat1[i, j];


				for (j = 0; j < i; j++)
				{
					mMat_r[i, j] = 0;
					for (k = 0; k < mLowerSize; k++)
						mMat_r[i, j] += (int)(((long)mTempMat1[i, k] * mMat_q[j, k]) >> PointPosition);
					for (k = 0; k < mLowerSize; k++)
						mMat_g[i, k] -= (int)(((long)mMat_q[j, k] * mMat_r[i, j]) >> PointPosition);
				}

				mMat_r[i, j] = 0;
				for (k = 0; k < mLowerSize; k++)
					mMat_r[i, j] += (int)(((long)mMat_g[i, k] * mMat_g[i, k]) >> PointPosition);
				mMat_r[i, j] = (int)(Math.Sqrt(mMat_r[i, j]) * Math.Sqrt(1<<PointPosition));

				for (k = 0; k < mLowerSize; k++)
					mMat_q[i, k] = (int)((long)mMat_g[i, k] << PointPosition / mMat_r[i, j]);
			}

			for (i = 0; i < mLowerSize; i++)
			{
				for (j = 0; j < i; j++)
				{
					mMat_g[i, j] = mMat_g[j, i] = 0;
					for (k = 0; k < i; k++)
						mMat_g[i, j] += (int)((long)-mMat_g[k, j] * mMat_r[i, k] / mMat_r[i, i]);
				}
				mMat_g[i, i] = 1<<PointPosition / mMat_r[i, i];
			}

			for (i = 0; i < mLowerSize; i++)
			{
				for (j = 0; j < mLowerSize; j++)
				{
					mMat_r[j, i] = 0;
					for (k = 0; k < mLowerSize; k++)
						mMat_r[j, i] += mMat_g[k, i] * mMat_q[k, j];
				}
			}

			#endregion

			// 更に右から転置行列を掛けます
			if (mRow < mCol)
			{
				for (i = 0; i < mCol; i++)
				{
					for (j = 0; j < mRow; j++)
					{
						mInverse[i, j] = 0;
						for (k = 0; k < mRow; k++)
							mInverse[i, j] += mMat_r[k, j] * mValue[k, i];
					}
				}
			}
			else if (mCol < mRow)
			{
				for (i = 0; i < mCol; i++)
				{
					for (j = 0; j < mRow; j++)
					{
						mInverse[i, j] = 0;
						for (k = 0; k < mCol; k++)
							mInverse[i, j] += mMat_r[i, k] * mValue[j, k];
					}
				}
			}
			else
				return new Matrix16(mMat_r);

			return new Matrix16(mInverse);
		}

		#endregion

		#region 演算子定義

		/// <summary>
		/// 指定された行列の正符号行列を定義します。
		/// </summary>
		/// <param basename="mat">正符号を与えられる行列クラスのインスタンスを指定します。</param>
		/// <returns>指定された行列をそのまま返します。</returns>
		public static Matrix16 operator +(Matrix16 mat)
		{ return (Matrix16)mat.Clone(); }

		/// <summary>
		/// 指定された行列の負符号行列を定義します。
		/// </summary>
		/// <param basename="mat">負符号を与えられる行列クラスのインスタンスを指定します。</param>
		/// <returns>指定された行列の正負を反転させた行列を返します。</returns>
		public static Matrix16 operator -(Matrix16 mat)
		{ return mat.Negative(); }

		public static Matrix16 Negate(Matrix16 mat)
		{
			return mat.Negative();
		}
		public static Matrix16 Plus(Matrix16 mat)
		{
			return (Matrix16)mat.Clone();
		}

		/// <summary>
		/// 行列クラスの加算演算子を定義します。
		/// この演算では指定された行列の値は変化せず、新しいインスタンスに結果が格納されて返されます。
		/// </summary>
		/// <param basename="mat1">加算される行列を指定します。被加算・加加算には演算上で実質的な意味の違いはありません。</param>
		/// <param basename="mat2">加算する行列を指定します。被加算・加加算には演算上で実質的な意味の違いはありません。</param>
		/// <returns>加算した結果の行列インスタンスを返します。</returns>
		public static Matrix16 operator +(Matrix16 mat1, Matrix16 mat2)
		{ return mat1.Add(mat2); }

		public static Matrix16 Add(Matrix16 augend, Matrix16 addend)
		{
			return augend.Add(addend);
		}
		/// <summary>
		/// 行列クラスの減算演算子を定義します。
		/// この演算では指定された行列の値は変化せず、新しいインスタンスに結果が格納されて返されます。
		/// </summary>
		/// <param basename="mat1">減算される行列を指定します。</param>
		/// <param basename="mat2">減算する行列を指定します。</param>
		/// <returns>被減算行列<paramref basename="mat1"/>から加減算行列<paramref basename="mat2"/>を減じた結果の行列インスタンスを返します。</returns>
		public static Matrix16 operator -(Matrix16 mat1, Matrix16 mat2)
		{ return mat1.Sub(mat2); }

		public static Matrix16 Subtract(Matrix16 minuend, Matrix16 subtrahend)
		{
			return minuend.Sub(subtrahend);
		}

		/// <summary>
		/// 行列クラスの乗算演算子を定義します。
		/// この演算では指定された行列の値は変化せず、新しいインスタンスに結果が格納されて返されます。
		/// </summary>
		/// <param basename="mat1">乗算される行列を指定します。被乗算・加乗算には演算上で実質的な意味の違いはありません。</param>
		/// <param basename="mat2">乗算する行列を指定します。被乗算・加乗算には演算上で実質的な意味の違いはありません。</param>
		/// <returns>乗算した結果の行列インスタンスを返します。</returns>
		public static Matrix16 operator *(Matrix16 mat1, Matrix16 mat2)
		{ return mat1.Mult(mat2); }

		public static Matrix16 Multiply(Matrix16 multiplicand, Matrix16 multiplier)
		{
			return multiplicand.Mult(multiplier);
		}

		/// <summary>
		/// 行列クラスの除算演算子を定義します。
		/// この演算では指定された行列の値は変化せず、新しいインスタンスに結果が格納されて返されます。
		/// </summary>
		/// <param basename="mat1">除算される行列を指定します。</param>
		/// <param basename="mat2">除算する行列を指定します。</param>
		/// <returns>被除算行列<paramref basename="mat1"/>から加除算行列<paramref basename="mat2"/>を除した結果の行列インスタンスを返します。</returns>
		public static Matrix16 operator /(Matrix16 mat1, Matrix16 mat2)
		{ return mat1.Mult(mat2.Inverse()); }

		public static Matrix16 Divide(Matrix16 divident, Matrix16 divisor)
		{
			return divident.Div(divisor);
		}

		#endregion

		#region objectメンバ

		/// <summary>
		/// 特定の型のハッシュ関数として機能します。
		/// </summary>
		/// <returns>ハッシュ値を返します。</returns>
		public override int GetHashCode()
		{
			return base.GetHashCode();
		}

		/// <summary>
		/// 行列クラスの等価を定義します。
		/// </summary>
		/// <param basename="obj">このインスタンスと等価であるかどうかを試験する行列インスタンス。</param>
		/// <returns>このインスタンスが指定されたインスタンスと等価であればtrue、それ以外はfalseを返します。</returns>
		public override bool Equals(object obj)
		{
			Matrix16 test = obj as Matrix16;
			if (test == null || test.Row != mRow || test.Col != mCol)
				return false;
			int i, j;
			for (i = 0; i < mRow; i++)
				for (j = 0; j < mCol; j++)
					if (mValue[i, j] != test.mValue[i, j])
						return false;
			return true;
		}

		/// <summary>
		/// 行列クラスの文字列表現を取得します。
		/// </summary>
		/// <returns>このインスタンスの行列の文字列表現を返します。</returns>
		public override string ToString()
		{
			string ret = "";
			int i, j;
			for (i = 0; i < mRow; i++)
			{
				ret += "| " + mValue[i, 0];
				for (j = 1; j < mCol; j++)
					ret += ", " + mValue[i, j];
				ret += " |\n";
			}
			return ret;
		}

		/// <summary>
		/// 行列クラスの文字列表現を取得します。
		/// </summary>
		/// <returns>このインスタンスの行列の文字列表現を返します。</returns>
		public string ToString(string format)
		{
			return ToString(format, CultureInfo.CurrentCulture);
		}

		/// <summary>
		/// 行列クラスの文字列表現を取得します。
		/// </summary>
		/// <returns>このインスタンスの行列の文字列表現を返します。</returns>
		public string ToString(string format, IFormatProvider provider)
		{
			string ret = "";
			int i, j;
			for (i = 0; i < mRow; i++)
			{
				ret += "| " + this[i, 0].ToString(format, provider);
				for (j = 1; j < mCol; j++)
					ret += ", " + this[i, j].ToString(format, provider);
				ret += " |\n";
			}
			return ret;
		}

		#endregion

		#region ICloneable メンバ

		/// <summary>
		/// このインスタンスのコピーのインスタンスを作成します。
		/// </summary>
		/// <returns>このインスタンスと等価な行列のインスタンスを返します。</returns>
		public object Clone()
		{
			return new Matrix16(mValue);
		}

		#endregion

		#region 特殊行列検査

		/// <summary>
		/// ゼロ行列検査。
		/// </summary>
		/// <returns>ゼロ行列であればtrue, それ以外の場合はfalseを返します。</returns>
		public bool IsZero()
		{
			int i, j;
			for (i = 0; i < mRow; i++)
				for (j = 0; j < mCol; j++)
					if (mValue[i, j] == 0.0)
						return false;
			return true;
		}

		#endregion
	}
}
